Integrative systematics and evolutionary history of Berylmys bowersi (Mammalia, Rodentia, Muridae)

Abstract The Bower's Berylmys (Berylmys bowersi) is one of the largest rodent species with a wide distribution range in southern China and the Indochinese Peninsula. The taxonomy and evolutionary history of the B. bowersi is still controversial and confusing. In this study, we used two mitochondrial (Cyt b and COI) and three nuclear (GHR, IRBP, and RAG1) genes to estimate the phylogeny, divergence times, and biogeographic history of B. bowersi. We also explored morphological variations among the specimens collected across China. Our phylogenetic analyses indicated that the traditional B. bowersi contains at least two species: B. bowersi and B. latouchei. Berylmys latouchei was considered a junior synonym of B. bowersi distributed in eastern China, which is confirmed to be distinguishable at specific level because of its larger size, relatively larger and whiter hind feet, and several cranial traits. The estimated split of B. bowersi and B. latouchei was at the early Pleistocene (ca. 2.00 Mya), which might be the outcome of the combined effects of climate change in the early Pleistocene and isolation by the Minjiang River. Our results highlight the Wuyi Mountains in northern Fujian, China, as a glacial refugia during the Pleistocene and call for more intensive surveys and systematic revisions of small mammals in eastern China.


| INTRODUC TI ON
are large rodents in the family Muridae that are widely distributed across Indochinese Peninsula, north-east India, and southern China (Denys et al., 2017). They are patchily distributed in various habitats, including primary and secondary forests, and plantations. Berylmys was originally described as a subgenus of Rattus by Ellerman (1947) with R. manipulus (Thomas, 1916) and R. berdmorei (Blyth, 1851) in it. Yet, Ellerman (1961) assigned R. bowersi (Anderson, 1878) to the subgenus Stenomys, and regarded mackenziei (Thomas, 1916) as a subspecies of R. bowersi. Musser and Newcomb (1983) reviewed the taxonomy of Berylmys and promoted it to full generic status, containing four species: B. manipulus, B. berdmorei, B. bowersi, and B. mackenziei; a treatment that was widely accepted subsequently (Agrawal, 2000;Denys et al., 2017;Musser & Carleton, 2005).
The Bower's Berylmys (Berylmys bowersi; Musser & Carleton, 2005) is the largest and most widely distributed species of the genus and has been recorded in southern China, Vietnam, Myanmar, northern Laos, north-eastern Thailand, India, and the Malay Peninsula (Denys et al., 2017). It is reported to live in burrows and is a habitat generalist inhabiting temperate and subtropical montane forests, disturbed forests, and abandoned jhum, cultivated fields, and moist deciduous and evergreen forests (Musser & Newcomb, 1983). It is a destructive agricultural pest in cereal and a reservoir or vector of numerous zoonotic diseases including scrub typhus, leptospirosis, and cryptosporidiosis (Boey et al., 2019;Chaisiri et al., 2017;Chen et al., 2019). A study on the evolutionary history and dispersal of B. bowersi is therefore important for public health and biodiversity conservation. Taxonomy of B. bowersi is still controversial due to limited research. It was described by Anderson (1879) based on specimens collected from the Kakyhen Hills of Husa (= Hotha), Dehong, Yunnan Province in China. Six other scientific names have been applied to samples of the species: including Mus latouchei Thomas, 1897;Mus ferreocanus Miller, 1900;Rattus kennethi Kloss, 1919;Rattus bowersii lactiiventer Kloss, 1919; Rattus wellsi Thomas, 1921;and Rattus bowersii totipes Tien, 1966; but all of them are currently treated as junior synonyms of B. bowersi (Denys et al., 2017;Wei et al., 2021).
Recently, we collected a specimen of B. bowersi in Qingliang Mountains, Anhui, eastern China. Our preliminary molecular analysis revealed that this specimen has a relatively high genetic distance with B. bowersi from Yunnan Province in southwestern China . Such a high genetic distance suggested underestimated diversity in the genus Berylmys. In the present study, we collected Berylmys specimens across China and sequenced two mitochondrial and three nuclear genes. By analyzing morphometric data and integrating newly generated sequences with the GenBank data, we assessed the taxonomy, phylogeny, and evolutionary history of B. bowersi.

| Sampling
Animals were handled complying with the animal care and use guidelines of the American Society of Mammologists (Sikes et al., 2016), and following the guidelines and regulations approved by the internal review board of Anhui Normal University (with no special approval number), and with the permissions of local government authorities. Specimens and tissues were deposited at Anhui Normal University, Anhui, China, and Kunming Institute of Zoology (KIZ), Chinese Academy of Sciences, Kunming, China. Specimens were identified based on their morphology and distributions following Thomas (1897), Musser and Newcomb (1983), and Smith & Xie (2009). We obtained 41 Berylmys specimens across China ( Figure 1). The specimen information and location data are listed in were amplified using the primer pairs provided in Table 2. The PCR products were purified and sequenced in both directions using the BigDye Terminator Cycle Kit v.3.1 (Invitrogen) on an ABI 3730xl sequencer (Applied Biosystems). Sequences of all genes were edited using SeqMan 7.1. Each gene was aligned using MEGA 11 (Tamura et al., 2021) and then checked manually. All new sequences were deposited in GenBank (Accession numbers OQ160975-OQ161004, OQ161020-OQ161049, OQ168972-OQ168387, Table 1). In addition, 32 corresponding sequences of 16 specimens were also downloaded from GenBank (Table 1). Homologous sequences of Mus pahari, Rattus exulans, Niviventer cremoriventer, Apodemus agrarius,
Divergence times were estimated using the concatenated mitochondrial-nuclear genes (mtDNA + nDNA) in BEAST v2.6.6 (Bouckaert et al., 2014). Because the missed genes will lead to inaccurate estimates of branch lengths and divergence times, we only used those samples with more than three sequences. Data blocks were defined based on genes and codon positions and evolutionary models or partition schemes were estimated based on the Bayesian Information Criterion (BIC). Three fossil calibration priors were used to provide for estimates of divergence times: (1) the most recent common ancestor (MRCA) among the "core" Murinae (12.5 Mya; 95% HPD: 11.31-13.99 Mya), with a lognormal distribution prior (mean: 12.5, standard deviation: 0.07, offset: 0), so the median age was at 12.5 Mya and the 95% CI was 11.1-14.0 Mya; (2) the major split occurred in the Indomalaya between Rattini and the remaining tribes of Murinae (11.2 Mya; 95% HPD: 9.91-12.67 Mya), also with a lognormal distribution prior (mean: 11.2, standard deviation: 0.08, offset: 0), so the median age was 11.2 Mya and the 95% CI was 9.8-12.7 Mya (Aghova et al., 2018); and (3) the oldest fossil of B. bowersi, which dated at 1.8 Mya, with an exponential distribution prior (offset = 1.8, M = 0.6 [1.8 × 0.333]), so the median age was 2.22 Mya and the 95% CI was 1.8-3.6 Mya (Zheng, 1993). Each Bayesian analysis was composed of a random starting tree, a lognormal relaxed molecular clock model, and a birth-death tree prior. Each analysis was run for 100 million generations, sampling every 10,000 generations. The first 10% of the samples were discarded as burn-in. Convergence was assessed using Tracer v1.7 (Rambaut et al., 2018).

| Species delimitation
The uncorrected p-distance of Cyt b gene was calculated in MEGA 11 to make pairwise comparisons of genetic differentiation within and between different phylogenetic lineages. BPP v3.1 was used for molecular species delimitation, employing both A10 (Yang & Rannala, 2010) and A11 models (Yang & Rannala, 2014). Based on the phylogenetic trees, we assigned B. bowersi and Berylmys sp. as two candidate species and used the best BI tree as the guide F I G U R E 1 Sampling localities of specimens used in the phylogenetic analysis.
TA B L E 1 Information of samples used for phylogeny analysis. in the Bayesian phylogenetics and phylogeography (BPP) analysis.

Species
To avoid poor mixing during BPP analyses, we followed the suggestions of Yang (2015) and used the combinations of algorithm 0, algorithm 1, and different priors for models A10 and A11, respectively. "Heredity = 1" allowing θ (referred to as "thetaprior") to vary among loci, or "Locus rate = 1" specifying the random-rates model of Burgess and Yang (Burgess & Yang, 2008), were also used but not at the same time. In total, 24 individual BPP analyses were performed (Table 4), and each was repeated two times.

| Morphological measurements and statistical analyses
External measurements including head and body length (HB), tail length (TL), hindfoot length without claws (HF), and length of ear from intertragal notch to crown (EL) were taken in the field with a ruler to the nearest 0.1 mm. The body weight (W) of each specimen was weighed to the nearest 0.01 g using an electronic scale.
Twenty-three craniodental measurements were taken to the near-  Thomas (1897), and Musser and Newcomb (1983) were followed for the terminology of morphological descriptions. All the craniodental measurements were taken by Zhongzheng Chen.
To evaluate the morphological variation between B. bowersi and Berylmys sp., variances of different variables between groups were tested using an independent sample t-test in SPSS 22.0, and a principal component analysis (PCA) was performed using the log 10transformed craniodental measurements.

| Phylogenetic relationships, divergence time, and species delimitation
We obtained 4950 bp long sequences for each specimen, including Berylmys sp. and B. bowersi on mitochondrial Cyt b gene sequence is 8.51% (Table 3). The BPP results consistently supported that Berylmys sp. and B. bowersi are two distinct species under different algorithms, models, and priors (PP = 1.00) ( Table 4).

| Morphological analyses
The external and skull measurements are given in Table 5. The mean values of most of the measurements of Berylmys sp. are larger than TA B L E 2 Primers and PCR conditions for amplification and sequencing used in the genetic analyses.  (Table 6). The first principal component (

| Morphological diagnosis
Because